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ABSTRACT 

Bridging the gap between first principles methods and empirical schemes, the density functional based 
tight-binding method (DFTB) has become a versatile tool in predictive atomistic simulations over 
the past years. One of the major restrictions of this method is the limitation to local or gradient 
corrected exchange-correlation functionals. This excludes the important class of hybrid or long-range 
corrected functionals, which are advantageous in thermochemistry, as well as in the computation of 
vibrational, photoelectron and optical spectra. The present work provides a detailed account of the 
implementation of DFTB for a long-range corrected functional in generalized Kohn-Sham theory. 
We apply the method to a set of organic molecules and compare ionization potentials and electron 
affinities with the original DFTB method and higher level theory. The new scheme cures the significant 
overpolarization in electric fields found for local DFTB, which parallels the functional dependence in 
first principles density functional theory (DFT). At the same time the computational savings with 
respect to full DFT calculations are not compromised as evidenced by numerical benchmark data. 

Keywords: Density functional based tight-binding ■ DFTB ■ long-range corrected exchange- 
correlation functionals • range separation ■ photoelectron spectroscopy ■ quasi-particle energies ■ 
static polarizability 


1. INTRODUCTION 

Density functional theory (DFT) is currently the de-facto standard in computational chemistry and 
condensed matter physics.® Due to algorithmic improvements and the availability of high performance 
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computing architectures, simulations with several thousands of atoms are becoming possibleP This 
development paves the way for in-silico combinatorial chemistry and materials scienceP® 


Notwithstanding, more approximate methods like the well established density functional based 
tight-binding method (DFTB) still have their meritsP The search for global minima on complex 
potential energy surfaces or the need for accurate phase space sampling in molecular dynamics 
simulations are just two examples in this contextP On the methodological level, DFTB has been 
extended in various ways during the last years. Such developments include the inclusion of spin-orbit 
interactions and non-collinear spin,l3 van-der-Waals correction^ and hybrid QM/MM schemes,l3E3 as 
well as the description of electronic excited states using the time-dependent formalism TD-DFTB.I^^^ 


Surprisingly, one of the most trivial steps in a DFT calculation, the change of the employed 
exchange-correlation (XC) functional, is comparatively cumbersome in DFTB. This is because DFTB 
relies on precomputed matrix elements of the Hamiltonian that are transformed to the molecular 
frame by Slater-Koster rules.l^ In addition, a consistent set of transferable pair potentials must 
be created by combined DFT/DFTB calculations on reference structures. While these steps are 
labor-intensive but straightforward for local and semi-local XC functionals, additional complications 
arise for hybrid and long-range corrected (LC) functionals that involve a fraction of Hartree-Fock 
exchange. In this case even the validity of the Slater-Koster rules is far from obvious and the quality 
of typical approximations in the DFTB framework like the neglect of three-center and crystal field 
terms must be newly assessed. 

In previous work, the formal foundations of DFTB with LC functionals were already outlined.^ 
In this contribution we provide details of an efficient implementation of these ideas and benchmark 
the accuracy of the new scheme. The interest in LC functionals arises due to their inherent reduction 
of the notorious self-interaction error in DFT. The resulting benefits of these functionals are well 
documentedThey yield the correct asymptotic -1/r behavior of the generalized Kohn-Sham 
potential which allows for the extraction of accurate ionization potentials and electron affinities (i.e., 
quasi-particle energies) from a single N-particle computation. In addition, the localization of the 
electron density in extended systems - the typical realm of DFTB applications - is strongly improved 
with respect to local functionals. We demonstrate that these desirable properties are also operational 
in the DFTB method which invokes additional approximations beyond the choice of the XC functional. 


After a specification of the XC functional used in this work in Section [27H we briefly summarize 
the results obtained in Ref. 16 before we provide a detailed account of our implementation. Section 
[^is devoted to the results of the new approach and contains a benchmark of fundamental gaps for a 
series of organic molecules. In addition we discuss the molecular response to electric fields (Section 


3.2) and comment shortly on the benefits of the new method for biological systems in Section 3.3 






2. THEORY 


2.1 Choice of the exchange-correlation functional 


In this work we use a functional that is based on the work of Baer, Neuhauser and Livshits.l^f^ 
Baer and Neuhauser use the adiabatic connection with a descreened electron-electron interaction 
v = [1 — exp {—ojr)]/r and derive an approximate form of the exchange-correlation functional. This 
form depends only on the range-separation parameter oj and provides the correct long-range behaviour 
of the potential: 
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The second term is given by 

N/2 

where Kohn-Sham molecular orbitals are denoted by xjji. The expression in Eq. approaches the 
exact Hartree-Fock (HE) exchange energy for large inter-electronic distance and is responsible for the 
reduction of the self-interaction error. 
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The explicit form of the DFT exchange-correlation functional has still to be approximated. 

Livshits and Baer suggested to use a correlation functional in the generalized-gradient approximation 
(GGA) together with short-range exchange in the local density approximation (LDA) according to 
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with some empirically determined parameter 0 < a < 1. The parameter set {a, a;} was subsequently 
optimized by fitting to thermochemical data. Similar empirical schemes are widely used to calibrate 
the whole plethora of different hybrid exchange-correlation functionals. l^l^ l^^HHI The aim of this 
work, however, is a proof of concept and we are interested in the qualitative improvement over the 
traditional DFTB scheme with respect to the failures due to the delocalization problem. Because of 
this we do not attempt any fine-tuning of the employed functional in the present work and simply 
use the PB^^Sl correlation functional and a = 1 in Eq.i 

For the LDA short-range DFT functional the exchange energy of the homogeneous 

electron gas with Yukawa interaction is usedP^HSU 
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Table 1: The compression radii [ag] for the elements used in this work. 


Results with this choice of xc functional will be termed LC-DFTB in the following. 


2.2 Total energy and Hamiltonian 

In this section we provide a brief overview of the LC-DFTB method, following Ref. [1^ together 
with an account of our practical implementation. The method is based on the linear combination of 
atomic orbitals (LCAO), where the single particle molecular orbitals -01 si'e expanded into a set of 
atom-centered basis functions 


Mr) = (7) 

As in the traditional DFTB method, a minimal valence-only basis set is obtained from pseudo-atomic 
DFT calculations with the additional confinement potential Konf(^") = {r/rM- The pseudo-atoms 
are considered to be spherically symmetric by equally distributing electrons over degenerate valence 
shells. In this study the functional is long-range corrected and defined by Eqs. [^to[^ The compression 
radius rg is usually taken to be proportional to the covalent radius of the respective atomic species. 
In the whole scheme, there are two confinement radii per element: amounts to the basis set 

compression and is used for the compression of the atomic densities in the initial guess.^^ For 

LC-DFTB we use the same confinement radii as previously reported for the mio-1-1 set,l^l^ which 
are summarized in Table We note that the basis set for sulfur features additional polarization 
functions.® 


The non-interacting reduced density matrix for the closed-shell case is given by 
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which defines the density matrix in the atomic orbital (AO) representation. 

With these notations and the particular choice of we are ready to write the total energy 

of the closed-shell LC-DFT in the AO basis 
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where single particle Hamiltonian hnv and two-electron repulsion integrals are given by 
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here denotes the external potential. Next, the density matrix is decomposed as P = Pq -l- AP, 
with some reference density matrix or initial guess Pq. Such a decomposition is common for direct 
SCF approaches where the Hamiltonian is constructed iteratively.^ We choose the reference density 
matrix to be the superposition of atomic density matrices Pq = which are obtained from 

the previous atomic LC-DFT calculations for each element. This choice parallels on the one hand 
the protocol in standard DFTB, on the other hand it is widely used as the initial guess in quantum 
chemistry codes. 

We now expand the local part of the exchange-correlation functional around the reference density 
matrix up to second order in AP: 
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First and second derivatives of the functional in the AO basis are denoted here as and 
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respectively. Eq.[l4| represents the first approximation in LC-DFTB, which amounts to linearization of 







the local exchange-correlation potential. Inserting this into the total energy expression and following 
the procedure in Ref. we rearrange the terms such that 
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where the last term, the repulsive energy ii^repP depends only on the reference density. The second 
approximation is to replace the repulsive energy by a sum of fast decaying pair potentials as it is 
the case in standard DFTB.^This is a reasonable approximation also in the present range-separated 
formalism as shown in Ref. [I^ 


The Hamiltonian evaluated at the reference density 
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is treated in the two-center approximation 
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where the density matrix of the dimer -|- is constructed with the density compression radius 
^density already mentioned, the basis functions 0^ stem from a LC-DFT calculation with the 
compression radius and are used to construct both the zeroth-order Hamiltonian matrix elements 
and the overlap matrix = / 0^(r)0i,(r) dv. Thanks to the two-center approximation, matrix 
elements for all possible geometries can be constructed from a small set of high-symmetry integrals 
according to Slater-Koster rules.l^ These parameters are thus tabulated as a function of inter-atomic 
distance Rab = iRyi — R-bI- 


Next we approximate the two-electron integrals using the Mulliken approximation 
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This allows for reduction of four-center integrals by a sum of two-center integrals 
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Figure 1: The long-range integral 7 ''^ for the carbon-nitrogen interaction as a function of inter-atomic 
distance for different values of the range-separation parameter 00 . The gray dashed line denotes the 
1/r limit. The full-range 7 ^*^ as defined in appendix of Ref. 33 is depicted with an orange dashed line. 


In the spirit of the original DFTB method,® these integrals are then parametrized as 
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where a G A, (3 E B, and is some atom specific decay constant, which has to be defined. In a 
similar fashion the full range two-electron integrals are simplified to 
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To obtain the final prescription for the evaluation of these integrals we follow the reasoning of the 
original DFTB method. For the off-site elements A ^ B, the contribution due to the xc-kernel 
is assumed to vanish. Thus the full-range 7 -integral reduces to a two-center Coulomb integral over 
simple spherically symmetric Slater type functions. For this case an analytical result is available,® 
which is also used in the present scheme. Exchange-correlation contributions are fully taken into 
account for the on-site case A = B, as shown below. To evaluate the long-range 7 -integrals, we 
extend the analytical formula of Elstner et al.® to the case of Yukawa interaction (see appendix[^for 
details). As example, we plot 7 ^ 1 x 1 for the carbon-nitrogen interaction as a function of inter-atomic 

















distance for some values of the range-separation parameter oj in Figurej^ Additionally, the full-range 
7 '^'^ evaluated with the analytical formula from Ref. 33 is also depicted (orange dashed line). 

For the on-site elements A = B, the integrals have the form (details are given in appendix 
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For each element, the Hubbard parameter [U = d‘^E/drA) is defined as the second derivative of the 
energy with respect to the occupation of the highest occupied atomic orbital. The decay constants 
ta may therefore be determined by the requirement that LC-DFT and LC-DFTB yield the same 
Hubbard value for each species. In the presence of a Fock term this relation is different from that in 
traditional DFTB [Ua = S/lbr^). The Hubbard parameter of an atom in LC-DFTB is instead given 
by 
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where / is the angular momentum of the highest atomic orbital (see appendix for details). Eval¬ 
uating Ua using first principles LC-DFT for a given range-separation parameter oo, Eq.j^is solved 
numerically for ta- 

Having defined all elements of the method we proceed by applying the variational principle to the 


total energy in Eq. 15 This yields the Hamiltonian 
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and the generalized Kohn-Sham equations 
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In contrast to the original DFTB scheme, the full density matrix needs to be evaluated self-consistently 
rather than just the Mulliken charges. 














2.3 Implementation notes 

The scheme has been implemented in the development version of the DFTB+ code.^^ With respect to 
the original DFTB method, two significant changes are necessary. First, the zeroth-order Hamiltonian 
(Eq. [It) ) now also includes a contribution due to screened Hartree-Fock exchange. The necessary 
adaption of the numerical integration routines is outlined in appendix)^ Since the and S matrix 
elements are precomputed there is no overhead involved in actual calculations. 


The second change is related to the presence of an exchange term in the Fock matrix, given by 
the second line in Eq. [^ This part is performance-critical, since a naive implementation would lead 


to scaling, despite of the applied Mulliken approximation (Eq. 18). 

We therefore resort to techniques which are widely used in direct SCF calculations.^^S The Hamil¬ 
tonian is constructed iteratively 


W(P") = + AP) = W(P"-i) + W(AP), 


(28) 


such that the Hamiltonian at the n—th iteration is the sum of the Hamiltonian at the previous 
iteration plus a correction. Here the density matrix at n—th iteration is labeled P"^ and AP denotes 
the difference between density matrices at different iterations of the self-consistent cycle, not the 


difference between P and Pq as in Eq. 14 


We then exploit the fact that the matrix elements iT(AP) get smaller upon approaching conver¬ 
gence. In order to make efficient use of integral prescreening techniques, the exchange matrix in Eq. 
26\ is rewritten in the following form 
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For each Hamiltonian sub-block (C,D) and atom pair (A,B), an estimate for the quantities QOg, € 
C, rr e H is given by 
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where Sab = ^{\Sai3\) and AP” = max(|AP”^|). Imposing the cutoff criterion 
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Figure 2: The average time per SCF cycle for the polyacene oligomers versus the number of basis 
functions for the threshold parameters Cthreshoid = 10“^®, 10“®, 10“®. The gray dashed line gives an 
extrapolation for ideal quadratic scaling. The inset shows the mean absolute error in the eigenvalues 
in Hartree for the cases of Cthreshoid = 10“® and Cthreshoid = 10“®. Calculations have been performed 
on a single core of an Intel Core-i7 CPU. The execution time was measured by the Linux time utility. 


we decide whether the contribution of the Flamiltonian matrix element is significant and avoid the 


evaluation if possible. The double sum on the right side of Eq. 32 is absorbed in Cthreshoid- 


We analyze the efficiency and accuracy of this prescreening approach by performing benchmark 
calculations for the polyacene series (C' 4 „+ 2 -Ff 2 n+ 4 ) up to n = 150. Figure depicts the scaling of 
the method with respect to the basis size. The gray dashed line shows the extrapolated quadratic 
function t{n) = cri^, with f(no) equal to the CPU time at the smallest oligomer size no = 5 and 

*^threshold 10 


For Cthreshoid = 10“^®, prescreening and exact evaluation according to Eq. M lead to identical total 
energies within machine precision. Even with this tight threshold criterion, the expected quadratic 
scaling is achieved. For the finite threshold values of Cthreshoid = 10“® and Cthreshoid = 10“® the scaling 
remains quadratic, however the prefactor is reduced by a factor of 2-3. In the inset the mean absolute 
error (MAE) of the eigenvalues with respect to the exact evaluation is shown. For Cthreshoid = 10“® 
this error does not exceed 10“® Ha. 


The overall scaling of the method is determined by the calculation step with steepest scaling. 
Since the Hamiltonian matrix construction scales quadratically in our scheme, the scaling for very large 
















systems should be determined by the 0{N^) behavior of the diagonalizer. However, for the tested 
systems the time spent for diagonalization is negligible compared to the Hamiltonian construction. 
Thus we expect cubic scaling only for much larger system sizes, especially if effective Divide &l 
Conquer diagonalization routines are available like in the current DFTB+ version. 

3. RESULTS 

After the discussion of the main approximations and computational efficiency of the LC-DFTB 
method, we benchmark its predictive power. At this point we focus on the electronic structure 
at fixed geometry to highlight the advantages with respect to the original DFTB method in the 
computation of electron removal and addition energies. This also represents a necessary first step 
in developing a time-dependent formalism for long-range corrected functionals in the spirit of the 
TD-DFTB scheme. 

3.1 Quasi-particle energies 

It is well known that KS eigenvalues from local DFT are poor estimates for quasi-particle ener- 
gjesHHIlZHlII Especially the eigenvalue of the highest occupied molecular orbital (HOMO), if inter¬ 
preted as electron removal energy (ionization potential) underestimates the experimental ionization 
potential (IP) by several eV. Electron affinities estimated from energy of the lowest unoccupied molec¬ 
ular orbital (LUMO) are likewise prone to large errors. As a consequence the HOMO-LUMO gap is 
much smaller than the experimental fundamental gap or the one obtained from the GW approxima- 
tion.l^E2lE3 This challenges the ability of the local DFT to provide useful single-particle picture of 
the physical systems. 

Two major problems of local DFT have been identified. On the one hand, the exponential 
asymptotic decay of the KS-potential (instead of the correct — 1/r behavior) leads to underbound 
electrons and a wrong description of the long-range interaction. On the other hand, the correct 
asymptotics alone is usually not sufficient for the correct prediction of quasi-particle energies.!^ The 
correct total energy of a system with fractional occupation N + 5, where N is an integer and 
0 < 5 < 1 is real, has been shown to exhibit a linear dependence on particle number between N 
and iV -I- 1. This is the necessary and sufficient condition to obtain the correct IP and fundamental 
gap. The linearity condition has been directly connected to the many-electron self-interaction error 
(MSIE) exhibited by local exchange-correlation functionals,^^ which in a real space picture manifests 
as a delocalization problem.^ 

The LC-DFT restores the correct asymptotics of the potential and shows remarkably close agree¬ 
ment with the linearity condition.l ^ l ^ l^l^ As a consequence, at least frontier orbitals in LC-DFT 







can be interpreted as electron removal energies. We expect the LC-DFTB method to show similar 
qualitative improvement compared to DFTB. 

To investigate this point, we choose a set of organic molecules for which experimental ion ener¬ 
getics data is available. Among others this set includes a selection of compounds that are relevant 
for photovoltaic applications and have been studied in a similar context with LC-DFT® and GW® 
methods. The structural formulas of these molecules can be found in the supporting material. 

All geometries have been optimized at the traditional DFTB level with the mio-1-1 sell^l^ and 
are used by default for both first principles and DFTB/LC-DFTB calculations. For the calculations 
in this work, we choose a value of cu = 0.3 ag ^ for the range-separation parameter unless stated 
otherwise. We found this value to give reasonable results for the prediction of ionization potentials. 
Akinaga and Ten-no showed, that for a range-separation of Yukawa type the optimal value of the 
range-separation parameter is usually higher than that for error function based functionals.^ In their 
work the optimization of the parameter was carried out on atomization energies of the G2-1 set with 
the cc-pVTZ basis. Flowever, it is so far not clear whether these parameters can be directly transferred 
to LC-DFTB. From our calculations, we found values ranging from u = 0.5 ag ^ to a; = 0.75 ag as 
suggested in Ref. [2^ to be too large for the prediction of accurate ionization potentials. In this case 
LC-DFTB tends to systematically overestimate IPs. 

We include the gradient-corrected PB^^ll functional, the hybrid functional BSLYF^zHlol gnd the 
long-range corrected BN L^^ functional in the comparison. In general a basis set of triple-zeta quality 
with polarization functions is employed, while we also use smaller sets for the BNL functional to 
estimate the basis set dependence of the results. All first principles DFT calculations have been 
performed with the NWCFIEM-6.3 package.!^ 

First we compare the negative of the FIOMO eigenvalue to the experimental ionization potential^ 
for the mentioned set of molecules. The deviation A = |IPexp| — I^homoI foi" LC-DFTB, DFTB, PBE, 
B3LYP and BNL is shown in the top part of Figure]^ It is found that LC-DFTB is quantitatively in 
better agreement with experiment than the standard DFTB and first principle approaches with the 
PBE and B3LYP functionals. The MAE for LC-DFTB is 0.50 eV compared to BNL/3-21G (0.67 
eV), BNL/cc-pVDZ (0.47 eV), and BNL/cc-pVTZ (0.29 eV). Since the deviations for the local and 
hybrid functionals are much larger [B3LYP/cc-pVTZ (2.04 eV), PBE/cc-pVTZ (2.87 eV), DFTB 
(2.50 eV)], we conclude that the assets of long-range corrected functionals are still visible at the 
approximated LC-DFTB level. 

For some molecules, most notably thiadiazole and methane, the deviation from both experiment 
and BNL can however be quite large. This seems to be an effect of the basis set employed in the 
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Figure 3: (Top) Deviation of the negative HOMO eigenvalue (—chomo) fi'om the experimental 
ionization potential for different levels of theory. (Bottom) Deviation of the HOMO-LUMO gap from 
the reference (BNL/cc-pVTZ). 


pyrazine 





















LC-DFTB method. As already mentioned in Sectionwe did not optimize the compression radii for 
the new functional and use the parameters of the mio-1-1 basis set without changes. Furthermore 
the minimal basis of the DFTB/LC-DFTB method might not provide enough variational flexibility. 
This can be deduced from the results for the BNL functional at the 3-21G, cc-pVDZ and cc-pVTZ 
level. Thiadiazole is an exception, since for this compound basis set effects in the first principles 
calculations (BNL/3-21G compared to BNL/cc-pVTZ) are small. Inspection of Figure]^ reveals that 
the results of the LC-DFTB method are in general comparable to LC-DFT results with small double- 
zeta basis (3-21G). The new scheme clearly outperforms first principles DFT calculations based on 
the PBE and B3LYP functionals for the description of ionization potentials. 

Next, we investigate fundamental band gaps. Since the experimental data for electron affinities 
is in general not available, we use BNL/cc-pVTZ results as reference. Deviations of the fundamental 
gaps from this reference are depicted in the bottom part of Figurej^ The MAE deviation for LC-DFTB 
is 1.36 eV, compared to DFTB (5.06 eV), BNL/3-21G (0.41 eV), BNL/cc-pVDZ (0.15 eV), PBE/cc- 
pVTZ (5.29 eV), and B3LYP/cc-pVTZ (3.92 eV). For the case of methane and dimethylether rather 
large deviations of 12.30 eV and 8.66 eV, respectively, are found for LC-DFTB. Again we assign these 
failures to the minimal basis set, since BNL/3-21G shows qualitatively similar, although much smaller, 
deviations of 3.58 eV and 1.57 eV for these molecules. In line with this we note that dimethylsulfide, 
which is essentially dimethylether with oxygen being replaced by sulfur, shows a much smaller error. 
This can be attributed to the fact that the sulfur in present parametrization contains additional 
polarization functions (d-orbitals). Inclusion of polarization functions for oxygen and nitrogen might 
thus reduce the mentioned problem with a moderate loss of computational efficiency. 

We now proceed by investigating the electronic structure beyond the frontier molecular orbitals 
which is experimentally accessible by photoemission spectroscopy. Here successful applications of 
theoretical methods such as GW or hybrid DFT have been reported, whereas local DFT exhibited 
serious flaws.l^l^^HH! The deficiencies of local DFT, which could be attributed to the self-interaction 
error, seem to be partially cured by long-range corrected functionals.l^®!^ The eigenvalue spectrum 
from standard DFTB, LC-DFTB, PBE/cc-pVTZ and BNL/cc-pVTZ theories for 3,4,9,10-perylene- 
tetracarboxylic-dianhydride (PTCDA) and pentacene (5-acene) molecules is presented in Figure]^ All 
spectra are rigidly shifted such that HOMO position is at 0 eV. We use a simple gaussian broadening 
profile with the full width at half minimum of 0.1 eV to mimic the experimental resolution and 
broadening. For this study geometries optimized at the B3LYP/cc-pVTZ level have been used. 

The experimental photoemission spectrum of the PTCDA molecule is characterized by the large 
gap of 1.5 eV between the first and second peaks, where the second peak appears at energies relative 
to the HOMO between -1.5 eV and -2.1 eV.^^ Long-range corrected functionals and GW are usually 
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Figure 4: Broadened density of states for PTCDA and pentacene molecules from LC-DFTB and 
DFTB compared to BNL/cc-pVTZ and PBE/cc-pVTZ results. The HOMO level has been shifted 
to the zero of energy for all methods. 


capable to reproduce these features qualitatively.^ We confirm this again with our BNL/cc-pVTZ 
calculation, see part a) of Figure]^ The spectrum is comparable to the standard LC-cuPBE result 
of Ref. ^ Since PBE provides poor results as has been already discussed in Refs. ^ 46 we do 
not expect a better performance for DFTB, which is essentially an approximate DFT with the PBE 
functional. As can be clearly seen from part c) and d) of Figure]^ the DFTB spectrum differs even 
more strongly from BNL/cc-pVTZ than PBE. 


We keep this fact in mind and proceed to the LC-DFTB results. We observe again a major 
difference to the BNL/cc-pVTZ spectrum. The LC-DFTB shows four cr-orbitals, mostly located at 
the anhydride groups of the PTCDA, right in the middle of the mentioned gap. These orbitals fit into 
the level ordering scheme of DFTB, therefore it seems that LC-DFTB just shifts the levels of DFTB 
if the range-dependent HF exchange term is added. However, we find for other systems that the 
shift is non-uniform and the level ordering is not preserved (see below the discussion of pentacene). 
For long-range corrected functionals the transition from a; = 0 (DFT limit) to finite uj values is 
quite generally accompanied by smooth inhomogenous shifts of the single particle levels. Therefore 
level reordering is expected. Obviously the BNL/cc-pVTZ theory correctly exhibits this reordering of 
levels with respect to PBE/cc-pVTZ, which is not observed for LC-DFTB in the example shown here. 
This should be attributed to the typical DFTB approximations, like the reduced basis set and the 












two-center approximation, since already the DFTB results in Figure]^) show strongly underbound 
cr-orbitals. As level ordering issues have also been observed in first principles LC-DFT, there can also 
be an additional effect. The analysis in Ref. 46 showed, that even if the frontier orbitals seem to be 
well described by LC-DFT, there are states, usually of different symmetry (e.g. cr-orbitals), which 
exhibit considerable orbital SIE (OMSIE). In the aforementioned work the spectrum of tuned LC- 
cjPBE functional has been discussed. Within the tuning procedure the value of the range-separation 
parameter is chosen such that the HOMO eigenvalue is equal to the ionization energy obtained from 
the difference of total energies of the neutral species and the cation. For this functional the second 
peak in the PTCDA spectrum is composed of degenerate a states (which correspond to the HOMO- 
2/HOMO-3 in the LC-DFTB spectrum in the present work). Analysis of the orbital self-interaction 
error for this theory showed large OMSIE in the cr-orbitals, while for 7r-orbitals the error was small. 
LC -cuPBE with standard value of the range-separation parameter exhibited the opposite behavior, 
where the cr-orbitals had rather small OMSIE. Thus the failure of LC-DFTB could also be connected 
to residual self-interaction error which is more pronounced for the anhydride cr-orbitals. In fact, 
the energetical position of the tt orbitals relative to the HOMO level is quite well represented by 
LC-DFTB, much better than in first principles PBE. 


The problematic level ordering of LC-DFTB is also observed for the pentacene spectrum (right 
part of Figure]^. While the levels up to HOMO-4 in DFTB, PBE and BNL show the same order, the 
LC-DFTB spectrum is characterized by the appearance of two cr-orbitals at HOMO-3 and HOMO-5 
positions. They are indicated by red and blue lines respectively. It turns out that these orbitals have 
already a differing position in the DFTB spectrum as compared to PBE. The approximate theories 
tend to underbind these orbitals. Thus in this case the influence of the DFTB approximations is 
more evident. From the analysis above, we conclude that the level ordering problem of the LC-DFTB 
might be partially caused by the orbital self-interaction error of the cr-states within the LC-DFTB 
theory as well as by the applied DFTB approximations and the minimal basis set used. Again we 
find a significant spreading of the quasi-particle spectrum in line with the LC-DFT results and quite 
accurate level positions for the HOMO, HOMO-1 and HOMO-2 levels of vr-nature. Notwithstanding, 
a full characterization of an experimental photoemission spectrum seems to be too ambitious at this 
point. 


At the end of this section we briefly comment on the computational efficiency of the new scheme 
versus first principles approaches. In fact, the main motivation to use approximate methods like LC- 
DFTB is the possibility to investigate large systems well beyond the scope of conventional DFT codes. 


In Section we already documented the quadratic scaling of the method with system size. In Table 
[^absolute timings of single point calculations for the theories considered in the previous discussion 




molecule 

BNL/3-21G 

BNL/cc-pVDZ 

BNL/cc-pVTZ 

LC-DFTBi 

LC-DFTB^ 

DFTB 

5-Acene 

325 
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11 

6 

1 

Perylene 

541 

4447 

113414 

12 

6 

1 

H 2 P 

932 

7486 

23825* 

16 

8 

1 

Coronene 

2116 

15983 

260677 

14 

6 

1 

6-Acene 

507 

3516 

79993 

15 

8 

2 

5-Thiophene 

1303 

11581 

144735 

16 

11 

2 

PTCDA 

3748 

20477 

521424 

22 

10 

2 

H 2 PC 

3034 

27231 

366838* 

39 

14 

2 

H 2 TPP 

5735 

43304 

744967* 

72 

22 

4 

n 

0 

7221 

65393 

655789* 

121 

23 
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Table 2: Wall time [s] of DFT calculation versus LC-DFTB for molecules with more than 30 atoms 
from the considered set. The DFTB results are given for comparison as well. The asterisk denotes 
the parallel jobs on 12 CPUs, LC-DFTB^ was performed with Cthreshoid = 10“^® and LC-DFTB^ with 

^threshold 10 


are summarized. The LC-DFTB calculations, with Cthreshoid = 10“^® (denoted by LC-DFTB^) and 
Cthreshoid = 10“® (denoted by LC-DFTB^) have been performed on a single core of an Intel Core-i7 
CPU as before. The execution time was measured by the Linux time utility, where the user time in 
seconds has been collected. For the DFT calculations the NWCHEM-6.3 code was used in serial and 
parallel versions. The serial version of NWCFIEM-6.3 was executed on Intel Xeon 2.8GHz machines, 
while the parallel Jobs were distributed over 12 CPUs on a cluster. Wall times were extracted from 
the NWCHEM-6.3 output files. 

As expected, the first-principles calculations are computationally more demanding. Even the 
calculation with small basis set at the BNL/3-21G theory level is at least 30 times slower than LC- 
DFTB^ and 50 times slower than LC-DFTB^ for smaller molecules. We note that the threshold 
parameter Cthreshoid = 10“® gives the eigenvalues with MAE errors below 10“® eV (see Figure]^, 
thus this choice can be considered as accurate for practical calculations. For larger systems the gap 
in computational time between the LC-DFTB and first-principles calculations increases due to the 
quadratic scaling of the LC-DFTB method. At the same time traditional DFTB is at least an order 
of magnitude faster than LC-DFTB with Cthreshoid = 10“®. 







3.2 Electric field response: DFTB vs. LC-DFTB 


A well known deficiency of local DFT is the exaggerated response to an applied electrostatic field. 
This behavior have been attributed to the lack of a necessary non-local response term in the exchange- 
correlation functional.lS^^ Thus all local and semilocal exchange-correlation functionals fail to pro¬ 
duce the correct induced field, which is opposed to the applied electric field. As a consequence, local 
DFT leads to a wrong density distribution, characterized by a too strong separation of the induced 
charge. This results in a drastic overestimation of static polarizabilities, which gets stronger with 
growing system size. Difficulties amplify for the hyperpolarizability and second hyperpolarizability.l^Sl 
This problem has also consequences for the application of DFT in the field of molecular electronics. 
The underestimated HOMO-LUMO gap, lack of a field-counteracting term and delocalization of the 
density lead altogether to a flawed description of transport properties, such as conductance.^ We 


showed already in the section 3.1 that the LC-DFTB provides essentially better description of the fun¬ 
damental gap, which suggests the reduction of the delocalization problem as compared to traditional 
DFTB. Recently, Sekino et al. provided evidence that LC-DFT shows the tendency to overcome the 
field response problem.ISS We therefore seek to confirm the signatures of a field-counteracting term 
due to the inclusion of the non-local range-dependent term in the LC-DFTB method. We calculate 
the static longitudinal polarizabilities of polyacetylene chains (PA, C 2 rM 2 n+ 2 ) with varying number 
of unit cells n and analyze the induced Mulliken charge distribution along the chain. LC-DFTB as 
well as traditional DFTB include the electric field F via the additional contribution 


-Enfield — — ^ Aq'aF ■ Ra (34) 

A 

to the total energy functional, where = Qa — Qa difference of the Mulliken population 

qa = J2u&a number of valence electrons on atom A located at position Ry^. We point 

out that traditional DFTB (and for the same reasons LC-DFTB) shows in general poor performance 
in predicting polarizabilities. This is due to the minimal basis set employed. However, LC-DFT even 
with minimal basis (e.g. BNL/ST0-3G) tends to correctly reduce the polarizability with respect to 
LDA/GGA-based theory in the same basis. 

We compare the polarizabilities obtained from the LC-DFTB method to long-range corrected DFT 
at the BNL/6-311G** and BNL/3-21G level. For the case of LC-DFT the polarizabilities have been 
obtained from coupled-perturbed Kohn-Sham theory (CPKS), implemented in the NWCHEM-6.3 
package. The geometries for PA with n = 10 and tt, = 40 have been optimized at the B3LYP/6- 
311G* level of theory. We obtain the LC-DFTB polarizabilities by applying the finite field method. 
The numerical derivative of the dipole moment /x along the long axis with respect to the perturbing 
electric field F is calculated with the center difference formula a = (/i(F) — /2F, where the 
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Table 3: Static longitudinal polarizability of PA (n = 10). All values are in atomic units [a^e^/ 


field strength was chosen to be F = 0.0004 a.u.. Tablelists the longitudinal static polarizability of 
the PA oligomer with n = 10 for BNL/6-311G**, BNL/3-21G and LC-DFTB at different values of 
the range-separation parameter cj. We observe similar qualitative behavior, although the quantitative 
differences are rather large, especially for the local DFT limit [u —)■ 0). In this local DFT limit all 
three theories exhibit larger polarizabilities as in the opposite HF-|-c limit (cu —)■ oo). In the latter case 
the exchange-correlation functional is composed of 100 % HF exchange and local DFT correlation. 
Remarkable is the rapid drop of the polarizability with the increase of the range-separation parameter 
oj (compare also Figure]^, seen in all theories. While in the case of n = 10 the ratio of HF-|-c to local 
DFT limits is 0.63 for both BNL/6-311G** and BNL/3-21G, and 0.68 for LC-DFTB, it decreases to 
0.23 for BNL/3-21G and 0.33 in the case of LC-DFTB for the larger system with n = 40 units. This 
indicates the aforementioned growing of the polarizability overestimation with the increasing system 
size. 

Further information is obtained by inspection of the charge density. The inset of Figure [^depicts 
the induced Mulliken charge due to a field of magnitude F = 0.001 a.u. along the n = 40 oligomer 
for different values of the range-separation parameter u. The LC-DFTB(a; —)■ 0) shows an almost 
linear charge distribution, which indicates an overly large polarization. We note that the DFTB 
is quantitatively very close to the LC-DFTB in the DFTB limit a; —)■ 0, although the exchange- 
correlation functional is slightly different. Thus we do not show the DFTB result for brevity. Increase 
of the parameter oj gives rise to an effective screening of the electric field, which leads to the correction 
of the polarizability towards more physical values. The LC-DFTB polarizability as a function of the 
range-separation parameter is provided in the main part of Figurefor comparison. We emphasize 
the qualitative nature of the LC-DFTB results presented here. Quantitatively correct polarizabilities 
require in general large basis sets. An alternative solution is the exploitation of empirical correction 
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Figure 5: LC-DFTB linear polarizability of the PA (n = 40) oligomer as a function of range-separation 
parameter oj. The induced Mulliken charge per unit (F = 0.001 a.u.) for different values of oj is 
shown in the inset. 

methods.l^l^ 

3.3 Proteins in gas phase and solution 

Experimental techniques that allow for the non-destructive extraction of proteins in combination with 
structural analysis make it possible to study intramolecular interactions in the absence of a partic¬ 
ular solvent.^ Electrospray ionization (ESI)f^ together with mass spectrometry and ion mobility 
measurements permit gas-phase structure determination. Likewise, recent progress in the develop¬ 
ment of X-ray free-electron lasers hold the promise to resolve protein structure in vacuum at atomic 
resol ution.!^!^ From the theoretical side, computations using solvent models in addition to gas 
phase calculations may finally provide a way to understand the protein folding mechanism in different 
environments. 

Efforts in this direction may address the still controversial question whether peptides adopt the 
zwitterionic form known from aqueous solution also in the gas phase.^^^® DFTB, like other DFT 
approaches based on local xc-functionals, shows difficulties in the description of the zwitterionic state 
where long-range charge-charge interactions play an important role. In a recent study Nishimoto et 
al.l^ found that the DFTB self-consistency cycle failed to converge for the model peptides chignolin 
(PDB ID: lUACPS) and Trp-cage (PDB ID: 1L2Y23I) in the zwitterionic conformation. It has been, 
however, shown both theoretically with force-field based molecular dynamics and experimentally using 

















Figure 6: The SCF convergence plot of chignolin and Trp-cage proteins as a function of range- 
separation parameter u for the LC-DFTB method. The SCF does not converge in the local DFTB 
limit {oj = 0). Structures of the proteins are shown as insets. These plots have been generated from 
PDB structuredly with the VMD software.l^ 

electrospray ionization and photo-dissociation, that the native zwitterionic configuration of Trp-cage 
remains stable in the gas-phase.l^l^ 

Raising the electronic temperature seems to circumvent the issue as a technical workaround, 
but in fact cures the symptoms instead of solving the fundamental problem, which is the notorious 
underestimation of the FIOMO-LUMO gap of traditional DFTB. Motivated by this we present single¬ 
point calculations on the lUAO and 1L2Y structures using the LC-DFTB method and compare to 
the LC-PBE/3-21G and RHF/3-21G methods in the following. The geometrical structure is in both 
cases the native zwitterionic conformation, obtained by NMR measurements as deposited in the 
protein data base.l^^ For our gas phase simulations we consider the models to be charge neutral in 
total. To this end basic and acidic side chains have been restored to their neutral form by appropriate 
protonation or deprotonationThe carboxyl and amino terminals were kept oppositely charged. A 
ribbon representation of both proteins is given in Figure]^ 

The LC-DFT and RFIF calculations have been performed with the parallel version of NWCFIEM- 
6.3 on 8 CPUs, whereas LC-DFTB was run on a single core. As can be seen from the Table a 

^The geometry of the affected functional groups was subsequently optimized at the DFTB level (electronic 
temperature T = 500 K) leaving the rest of the molecule fixed. 
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LC-DFTB(a; ^ oo) 


78 
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8841 
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Table 4: Frontier orbital energies and fundamental gap (all in eV) of the chignolin and Trp-cage 
zwitterions in the gas-phase for different theories. 


good agreement in the description of the frontier orbital energies between full DFT and LC-DFTB is 
achieved. The timings can be found in the bottom part of the Table. We also performed LC-DFTB 
calculations for different values of the range-separation parameter oj. The convergence (number of 
SCF cycles) for both proteins is depicted in the figure]^. As can be clearly seen, the convergence 
issue does not occur for typical values of the range-separation parameter. As the band gap opens, 
convergence generally improves. 

These results indicate that LC-DFTB allows for reliable studies of proteins in the gas-phase both 
in the neutral and zwitterionic conformation. Certainly, further studies need to address also the 
energetics of the various conformers but already at this point LC-DFTB seems to be a useful tool for 
biological applications. In this context it might serve as the underlying electronic structure method 
for the fragment molecular orbital approach (FMO),®® which enables the quantum chemical study 
of biological systems with many thousands of atoms. 

4. CONCLUSION AND OUTLOOK 

In the present paper we addressed the implementation and benchmark of a long-range corrected 
exchange-correlation functional in the DFTB method, which we denote as LC-DFTB, for closed- 
shell systems. The practical implementation requires the extension of tools for the evaluation of 
Slater-Koster parameters as well as some rather minor changes to the DFTB Hamiltonian. Thus the 
scheme can be easily implemented in existing ab Initio software packages, which in principle contain 
all necessary routines for the integral evaluation and Hamiltonian construction. 










The Hamiltonian construction exhibits quadratic scaling and dominates the computational time 
for the systems tested. For larger systems the method is expected to show cubic scaling due to 
the diagonalization. The performance benefit with respect to first principles DFT with long-range 
corrected functionals still remains two to three orders of magnitude depending on the size of the 
basis set. 

LC-DFTB shows clear signatures of correction of the delocalization problem in local DFT that 
is attributed to the self-interaction error. This has been demonstrated for the frontier orbitals of 
a set of organic molecules, where in general a promising agreement with full LC-DFT with double 
zeta basis (BNL/3-21G) has been observed. Remaining flaws were related to the minimal basis set 
characteristic for DFTB, which influences especially electron affinity levels. As a second example we 
calculated the static longitudinal polarizability of polyacetylene chains and provided evidence for the 
qualitatively correct description of the response potential. 

The parametrization of the repulsive potential, which requires an adjustment due to the modified 
Hamiltonian, is presently on the way. With this development, thermochemical and structural prop¬ 
erties may be investigated in addition to the electronic structure. The extension of the method to 
a spin-unrestricted formalism is likewise promising, since it would allow for the non-empirical tuning 
of the range-separation parameter.^ Finally we would like to point out that the presented formalism 
is not restricted to the specific long-range corrected functional used in this study. The scheme can 
be easily adapted to conventional hybrid functionals like B3LYP or to screened exchange functionals 
like HSE,E! which are recently becoming increasingly popular for periodic systems and solid state 
applications. 
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APPENDIX A. NOTE ON THE NUMERICAL EVALUATION OF INTEGRALS 

The evaluation of electron repulsion integrals for Slater type orbitals (STO) as used in the DFTB 
method is a numerically challenging task. For the Coulomb interaction several algorithms based on 
the analytical expansions of STOs around a center are available. Recently, Seth and Zieglei® ex¬ 
tended these techniques to the case of a Yukawa interaction (exp(—ci;r)/r). In the present method 
the run-time evaluation of integrals is explicitly avoided by invoking the Slater-Koster scheme with 


precomputed parameters. Thus the integral evaluation does not affect the computational perfor¬ 
mance. We therefore resort to a simple, yet robust, numerical integration scheme proposed by 
Becke.®^ This scheme was designed for the evaluation of two-electron integrals over numerically 
defined charge distributions with Coulomb interaction. The modification for the case of Yukawa 
interaction is straightforward and laid out below. 

The double integral over four arbitrary orbitals 0j(r), which are in general located at different 
centers, and some interaction g{\r — r'|) 

{ab\cd) = j j 0a(r)06(r)5'(|r - r'|)0c(r')0d(r') dvdr', (35) 


can be seen as a special case of the integral over general charge distributions p(r), cT(r) 


/ = 



p{r)a{r')g{\r -r'\) drdr' = / p(r) 


^(r')5'(|r - r'\)dr' 


dr. 


(36) 


VM 

The potential fC(r) can be obtained by solving the following boundary value problem with fC(r) —>■ 
0, |r| —> oo 


Wir) = — 47rcr(r) 


(37) 


with some differential operator V. The function p(|r — r'|) 
problem 


in eq. 


36 is the Greens function of the 


'Dg{\r — r'l) = —47r(5(r — r'). (38) 

The differential operator for the Coulomb interaction g{r) = l/r \sV = V^. The Yukawa interaction 
g{r) = exp(— a;r)/r is the Greens function for the modified Helmholtz operator P = — o;^. Thus 

the Poisson equation in the original work by Becke has to be replaced by the modified Helmholtz 
equation. 

To make use of efficient quadrature schemes, Becke proposed a space partitioning and decompo¬ 
sition of the charge distribution into atom-centered portions. The space is partitioned by some set 
of functions /A(r), with = 1, Vr G Using this set of functions the density is divided 

into the atomic contributions 


pw = W- (39) 

A A 

^The choice of the numerical integrator can be additionally motivated by the fact that the DFTB method 
does in principle not rely on a specific type of orbitals as basis functions. This allows for more flexibility, 
especially in the context of atomic DFT calculations with converged numerical orbitals. 








This allows to write the integral in Eq. 36 as a sum of two-center integrals over atom-centered charge 
distributions pA, cb 


/ = ^ /jB = ^ / / PA(r)<rB{i')9{\i - r'l) **'. 


(40) 


AB AB 

The inner part of the two-center integral Iab is evaluated as follows. In the Helmholtz equation 


(V^ Va{y) = -AnpAir), 

the potential and density are expanded into spherical harmonics 

va{v) = 

Im 

PaM = ^PAJm(r)yim((l)- 


(41) 

(42) 

(43) 


Im 


Inserting these expansions into Eq. 41 gives a set of one-dimensional radial equations 


^ _ /(/ + !) 

dr‘^ 


UJ 


Vb,;■»(>■) = -4irr/3B,im(>'), 


(44) 


with boundary conditions Yimr^QVA,im{f') = 0, \im.r^ooVA,im{f') = 0. These equations are solved 
using a finite difference method and from the resulting VA,im the potential VA(r) is assembled 


according to Eq. 42 The remainder of the integration is performed on a two-center grid as described 


in the original paper by Becke.^ 

APPENDIX B. ANALYTICAL EVALUATION OF THE LONG-RANGE 7 -INTEGRAL 

The 7 -integral over the Yukawa interaction can be reduced to a one-dimensional integratiorP^ 
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where Rab = IR-a — Rb|- This integral can now be further evaluated using the residue theorem. 
We obtain then the analytical expressiorj^ 
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In the limit a; —)■ 0 the integral is taken over the Coulomb interaction. In this case our analytical 
formula reduces to the result of Elstner et al.^ The long-range 7-integral is then given by the 
difference i^B = Iab “ Iab- 

To obtain the on-site value [Rab —^0, A = B), we go back to the one-dimensional integral 
Eq. 46 We expand the term sin(gRAA) around zero and integrate again using the residue theorem 
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This yields the on-site value for the integral Eq. 
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APPENDIX C. THE HUBBARD PARAMETER IN LC-DFTB 

We derive the expression for the Hubbard parameter of atomic LC-DFTB. Note that in this case 
tU'' = 'ifA and El, = Using the orhogonality condition '^yCy^Cy^j = dij we 

obtain the total LC-DFTB energy of one atom in terms of occupation numbers n* 
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(52) 

(53) 


The analytical formula Eq. 47 has been derived after the calculations for this work have been done. 


For the presented results the long-range integral has been evaluated with a numerical integrator. We note 
that the analytical formula is more practicable and is recommended for use in practical implementation of 
the method. 
























We consider further only terms, quadratic in the occupation numbers, which belong to the highest 
occupied shell. It contains n electrons, equally distributed over the shell. The occupation number 
for an orbital is then n/di, where di = 21 + l\s the degeneracy of the shell with angular momentum 

1. Then the energy can be written as 
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(54) 
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Thus the second derivative with respect to the occupation numbers of the highest occupied shell is 

(57) 
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And together with Eq. 51 the atomic Hubbard from the LC-DFTB is obtained 
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This equation defines the decay constant r. 
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